Occupation numbers of the harmonically trapped few-boson system 
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We consider a harmonically trapped dilute A''-boson system described by a low-energy Hamil- 
tonian with pairwise interactions. We determine the condensate fraction, defined in terms of the 
largest occupation number, of the weakly-interacting A'^-boson system (A'^ > 2) by employing a per- 
turbative treatment within the framework of second quantization. The one-body density matrix and 
the corresponding occupation numbers are compared with those obtained by solving the two-body 
problem with zero-range interactions exactly. Our expressions are also compared with high precision 
ab initio calculations for Bose gases with N = 2 — 4 that interact through finite-range two-body 
model potentials. Non-universal corrections are identified to enter at subleading order, confirming 
that different low-energy Hamiltonians, constructed to yield the same energy, may yield different 
occupation numbers. Lastly, we consider the strongly-interacting three-boson system under spher- 
ically symmetric harmonic confinement and determine its occupation numbers as a function of the 
three-body "Efimov parameter" . 
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I. INTRODUCTION 



The weakly-interacting homogeneous Bose gas has 
been studied extensively in the literature 1-6]. Most 
commonly, the equation of state of the homogeneous Bose 
gas is expressed in terms of the square root of the dimen- 
sionless gas parameter p[as(0)]'^, where p denotes the den- 
sity and as(0) the zero-energy s-wave scattering length. 
The leading order term is the mean-field energy and the 
lowest order correction accounts for quantum fluctua- 
tions. In an alternative approach d, 0-l3i the ground 
state energy of A^ bosons in a cubic box with periodic 
boundary conditions has been obtained by applying per- 
turbation theory or the techniques of effective field the- 
ory to the weakly-interacting regime. As outlined by Lee, 
Huang and Yang [S] , the latter approach must reproduce 
the equation of state of the weakly-interacting homoge- 
neous Bose gas if the energies of the "subclusters" are 
summed up carefully. 

In addition to the energy, other observables of the ho- 
mogeneous weakly-interacting Bose gas have been con- 
sidered. The condensate fraction Nq/N, i.e., the frac- 
tion of particles in the macroscopically occupied low- 
est momentum state, is a particularly interesting quan- 
tity since it can be measured experimentally. Further- 
more, the connection between the condensate fraction 
and the superfluid fraction has been investigated in the 
literature, starting with the seminal works of London, 
Penrose and Onsager, and others |10l - |l3| . The conden- 
sate fraction has, as the energy per particle, been ex- 
panded in terms of the gas parameter p[as(0)]'^, Nq/N = 
1 — 8/(3^/7?) •\/p[as(0)]'^ + • • • . Application of the local 
density approximation shows that the condensate frac- 
tion of the weakly-interacting Bose gas under spheri- 
cally symmetric harmonic confinement scales as N/Nq — 
1 - 5v/7r/8v/p(0)[as(0)]3 + ■■■ , where p(0) denotes the 
peak density |14l |. 

This work considers A^ identical mass rUn bosons un- 



der spherically symmetric harmonic confinement with 
angular trapping frequency uj. For a review article of 
trapped gases, the reader is referred to Ref. [l^- In 
the weakly- interacting regime, i.e., in the regime where 
the two-body s-wave scattering length as{0) (expressed 
in units of Oho) and the product of the two-body effec- 
tive range Ve and [as(0)]^ (expressed in units of aj^^) 
are small, we determine expressions for the condensate 
fraction Nq/N; here, Oho denotes the harmonic oscillator 
length. Oho = ^/fi/i'tnaUj). For trapped systems, the con- 
densate fraction is related to the largest eigen value of the 
one-body density matrix. In particular, the largest eigen 
value or occupation number of the one-body density ma- 
trix defines the condensate fraction. Our results are ob- 
tained by applying time- independent perturbation theory 
to the A^-boson Hamiltonian with pairwise zero-range in- 
teractions characterized by 0^(0) and re[as(0)]^. The per- 
turbative expressions are compared with highly accurate 
numerical results for Bose gases with N ~ 2 — A that in- 
teract through a sum of short-range two-body model po- 
tentials. This comparison confirms that the leading-order 
term of the condensate depletion scales as (A^— l)[as(0)]^. 
At sub-leading order, a non-universal correction appears, 
i.e., a correction which is independent of as{0) and r^, 
and which is not needed to reproduce the energy of the 
finite-range system within an effective field theory ap- 
proach p^.[T7|. 

For the two- and three-boson systems, we go beyond 
the weakly-interacting regime. For two harmonically 
trapped bosons that interact through a regularized zero- 
range interaction potential, we determine the occupa- 
tion numbers as a function of the scattering length. For 
the three-boson system, we consider the unitary regime 
[l/as(0) = Te = 0] and determine the occupation num- 
bers as a function of the three-body phase or Efimov 
parameter. The occupation numbers for the two- and 
three-body systems show "oscillatory behavior" in the 
positive energy regime if plotted as a function of the rel- 



ative two-body energy and relative three-body energy, 
respectively. In the two-particle case, the oscillations are 
associated with the fact that the two-body s-wave phase 
shift changes by 27r as the two-body energy changes by 
about 2?kj. In the three-particle case, in contrast, the os- 
cillations are associated with the fact that the three-body 
Efimov phase goes through cycles of 27r as the three-body 
energy changes. 

The remainder of this paper is organized as follows. 
Sectionllllintroduces the system Hamiltonian and defines 
the one-body density matrix and the occupation num- 
bers. Section IIIII discusses the occupation numbers of 
the trapped two-boson system. Section |IV] considers the 
weakly-interacting regime of the iV-boson system. Sec- 
tion|V]considers the strongly-interacting three-boson sys- 
tem. Lastly. Sec. lVII concludes. Mathematical details are 
relegated to Appendix [^ and Appendix [Bj 



II. SYSTEM HAMILTONIAN AND 
DEFINITIONS 



scattering length is defined by taking the scattering en- 
ergy to zero, i.e., as(0) = limfe^oOsC^)- In many cases, 
the energy-dependence of as{k) is neglegible and as{k) 
can be replaced by the zero-energy scattering length 
as(0). In other cases (see Sees. IIIII and ITVT) . it is con- 
venient to parameterize the energy dependence of as(fc) 
in terms of the effective range r^, and the shape or volume 
parameter V [20], 



1 



^ irefc2 + ly/fc4^0(fc6) (5) 



a, (ft) a,(0) 2 



a,(fc) = a,(0) -t- ha,{0)]\,e ~ ha,{0)]''Vk^ + 0(fc«).(6) 

Z <5 

We note that as(0), r^ and V are only defined if the 
two-body potential falls off faster than r~^ , r~^ and r~J, 
respectively, in the large rj-fe limit i21n22i] . The pseudopo- 
tential given in Eq. ([3]) can alternatively be parametrized 
through the boundary condition [23l . |24| | 



We consider N identical mass nia bosons that inter- 
act through a short-range interaction potential Vtb un- 
der external spherically symmetric harmonic confinement 
with angular trapping frequency u. For this system, the 
Hamiltonian H reads 



N N 

J? = ^ffho(r,) + ^l^tb(r,fe), 



(1) 



where Hho{Yj) denotes the single-particle harmonic oscil- 
lator Hamiltonian, 



^ho(rj) 



2m, 



-^r, + 7:"^aW^r?, 



(2) 



and Vj the position vector of the j boson measured 
with respect to the center of the trap. We consider three 
different short-range model potentials Vtb(rjfe)7 where 



r.7fc 



r/c- 



Our two- and three-boson studies discussed in Sees. IIIII 
and |V] employ the regularized pseudopotential Vps [j| , 



yA^jk) 



^rfo^^O). ,^ 



dr 



'jk 



-rjk, 



(3) 



where rjk = \Yjk\- The operator {d/drjk)rjk ensures that 
the A'^-particle wave function tp is well behaved when the 
interparticle distance rjk goes to zero. In Eq. Q, as(k) 
denotes the energy-dependent scattering length ^8|, |l3] , 



as{k) 



ta,n{So{k)) 
k ' 



(4) 



where k denotes the wave vector associated with the 
scattering energy El°^ in the relative coordinate, k = 
y^rriaEl^^/h, and 6o{k) the energy-dependent s-wave 
scattering phase shift. The "usual" (zero-energy) s-wave 



3(ri2V'(ri2,Ri2,r3,--- .rjv)) 
dri2 

ri2'0(i"i2,Ri2,r3,--- ,rjv) 



ri2^0 



sikV 



(7) 



where R12 — (ri + r2)/2. The limit ri2 — > on the left 
hand side of Eq. (O is taken while keeping the coordi- 
nates R12, ra, • • • , Fat fixed. Analogous expressions hold 
for the other interparticle distances. 

In our perturbative calculations (see Sec. IIV[) . in con- 
trast, we write Vtb as a sum of the unregularized or bare 
Fermi pseudopotential Vp [25l |. 



VF(rjfe) = d'^-'>{rjk), 

and the zero-range potential V [tI, [l7J. 



(8) 



V'{v,k) = - 



Tih^[asm\, 



ylj^'\v,k)+5^^\v,k) 



which accounts for the effective range dependence. The 
first and second Laplacian in Eq. ([9]) act to the left and 
right, respectively, and ensure that the pseudopotential 
V' is Hermitian. While a pseudopotential that depends 
on the shape parameter could be added, it is not con- 
sidered here since it leads to higher order contributions 
in aj^" than we are considering in Sec. IIVI Since Vp and 
V' are not regularized, their use within perturbation the- 
ory leads to divergencies, which can be cured within the 
framework of renormalized perturbation theory by in- 
troducing appropriate counterterms denoted by W . We 
treat Vp and V' in second- and first-order perturbation 
theory (see Appendix 1X1). This implies that W must con- 
tain a term proportional to [as(0)]^ that cures the diver- 
gencies arising from Vp] no divergencies in the energy 
arise when treating V' in first-order perturbation the- 
ory[il,[i3. 
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FIG. 1: (Color online) Scattering quantities for Gaussian 
model potential Ig in the "weakly-interacting" and "strongly- 
interacting" regimes, (a) The solid line shows the quantity 
[as(0)]^re/ro as a function of aa(0)/ro (weakly-interacting 
regime), (b) The dashed line shows the quantity re/ro as 
a function of [as(0)/ro]~^ (strongly-interacting regime). Dot- 
ted lines are shown to enhance the readability of the graphs. 



Lastly, our numerical stochastic variational calcula- 
tions (see Sec. lIVp employ a short-range Gaussian po- 
tential Vg with range r^ and depth Vq, 



^g(^jfc) = Voexp 



V2ro 



(10) 



For a fixed rg, Vo is adjusted so as to generate potentials 
with different two-body s-wave scattering lengths as(0). 
Throughout, we limit ourselves to parameter combina- 
tions such that Vg supports no two-body bound states 
in free space. This implies that Vg is purely repulsive, 
i.e., Vo > 0, for 0^(0) > 0. For as(0) < 0, we have 
Vo < 0. The leading order and sub-leading order energy- 
dependence oi as{k), parameterized by fg and V, respec- 
tively, depend on rp. 

The solid line in Fig. [TJa) shows the quantity 
[as{0)]'^re/rQ [see Eq. ([6|)] for the Gaussian model po- 
tential as a function of the s-wave scattering length 
as{0)/ro. It can be seen that [as(0)]^re/rQ goes to zero 
as |as(0)/ro| — > 0. Moreover, [as(0)]^re/rg is positive 
for small negative as{0)/rQ and negative for small posi- 
tive as{0)/rQ. When |a5(0)| becomes infinitely large, the 
deviations from universality depend on r^ [see Eq. ([S|)]. 
The dashed line in Fig. [IJb) shows the effective range 



re/ro for the Gaussian model potential as a function of 
the inverse s-wave scattering length [as{0)/ro]~^- The 
effective range r^/rQ is finite and positive as |as(0)| — ^ oo 
and varies approximately linearly for small [as(0)/ro]^^ 
with negative slope. 

The higher-order energy-dependence of the s-wave 
scattering length in the "weakly-interacting" and 
"strongly-interacting" regimes is governed by the volume 
parameter V. The quantity [as(0)]^F/r3 [see Eq. ([6])] 
goes to zero as |as(0)| — >■ 0, and is positive for small 
as{0)/ro < and negative for small as{0)/ro > 0. The 
quantity V/r^ [see Eq. dS])] is finite and negative when 
the s-wave scattering length diverges; V/r^ varies ap- 
proximately linearly for small [as(0)/ro]^^ with positive 
slope. 

Sections IIIHTVl present results for the occupation num- 
bers Ui,, which are — for inhomogeneous systems — defined 
by way of the one-body density matrix p{r[,ri) [ll|, IZo - 



p(r;,ri)=iVx 

/[V^(r;,r2,--- ,rAr)]*i/;(ri,r2, 



(11) 



im 



,rjv)|2d3ri---d3rjv 



The one-body density matrix p(r']^,ri) can be expanded 
in terms of a complete orthonormal set. 



P(r'i,ri) 



E 



nt.0*(ri)0j,(r'i). 



[12) 



where u collectively denotes the three quantum numbers 
z/A/i needed to uniquely label the functions (/)^(ri) of the 
complete orthonormal set. If we use spherical coordi- 
nates, ly is the radial label, A the partial wave label and 
/i the corresponding projection number. The quantities 
0iy(ri) and n^, are called natural orbitals and occupa- 
tion numbers, respectively. Our normalization is chosen 
such that ^^ Tin = N . The largest occupation number 
rii, defines the condensate fraction Nq/N of the TV-boson 
system, i.e., Nq/N ~ iRax{n^/N). 

In practice, it is convenient to define partial wave pro- 
jections pxf_,{r[,ri), 

PXf.{r'i,n) ^ J Y^^{f[)p{r[,r,)Yx^{h)d'r[d'h, (13) 

where d^fi and d^f[ denote angular volume elements. 
The two-dimensional projections pxf_i (r'l , ri ) can be diag- 
onalized, yielding the occupation numbers n^, = Uvx^i, 
and the radial parts Pv\{ri) of the natural orbitals 
0iy(ri), where P„x{ri) is defined through 01/ (ri) = 
P.x{ri)Yx^,{h)- 



III. TRAPPED TWO-BODY SYSTEM 

The eigen energies and eigen states of the two-particle 
Hamiltonian are most easily determined by transform- 
ing the Schrodinger equation to center of mass and rel- 
ative coordinates R12 and ri2. In these coordinates. 



the wave function ip separates into the center of mass 
wave function i/jq^m (-^12) ^'^d the relative wave func- 



tion ipqfm{TCi2)- The two-body energy then reads 



where 



and 



E2 - ^r + ^2' 



E^"" = {2Q + L + 3/2)hu 



E!^"^ = {2q + l + 3/2)nu 



(14) 



(15) 



(16) 



with Q = 0, 1, • • • , L = 0, 1, • • • and ? = 0, 1, • • • . For 
each L (l), the center of mass (relative) energy has a 
2L + 1 {21 + 1) degeneracy that is associated with the 
projection quantum number M (m). The allowed values 
of q, and consequently the radial parts of the relative 
wave function and the relative eigen energies, depend on 
the functional form of the interaction potential Vtb- 

For the energy-dependent zero-range potential V^s [see 
Eq. ([3])], the relative wave functions with / > are not af- 
fected by the interaction potential, implying g = 0, 1, • • • ; 
in this case, the relative wave function coincides with the 
harmonic oscillator wave function and the two-body en- 
ergy is independent of the s-wave scattering length. For 
1 = 0, the relative eigen energies are obtained by solving 
the transcendental equation 29] 



V2r(3/4 - 



Ef/{2hw)) _ flho 



r(l/4-£;fV(2;ia;)) 



a,(ii;fi) 



(17) 



where the energy-dependent s-wave scattering length is 
evaluated at the relative energy of the trapped system, 
i.e., where we have set E^f = Sf' [Slil. Sohd hues in 
Fig.[2]show the relative eigen energies with I = obtained 
by solving Eq. pT)) with as(i?2°') replaced by as(0) as a 
function of as(0)/aho- In general, the eigen energies of 
the trapped two-body system need to be determined self- 
consistently since £^n°' ap pears on the left and right hand 
sides of Eq. (d?]) [ll[i3. 

For |a5(£'f')/aho| < 1, we Taylor expand Eq. p7)) 
around the non- interacting relative energies i?2 «; where 
^2% = i2n + 3/2)huj with n = 0, 1,---. Replacing 
l/as(i?2°') by the right hand side of Eq. dS]), we find 



E. 



rcl 



V- ,{^.]) i as{0) \ I r, 

^2,n T" /_^ "2,n 

i=i,j=0,j<i 



Ef_. + 



flho 



flho 



f)W 



.(18) 



The next terms are proportional to [as(0)]^, [as(0)]*re, 
[as(0)]3r2 and [as{Q)fV. The coefScients S^'.n can be 
compactly written in terms of the quantity hn^p, 



hn,p — Hnp + (—1) -ff-n-3/2,pi 



(19) 



where Hn,p is a generalized harmonic number [30|. Ex- 
plicit expressions for the coefficients d^n with i + j < 4: 






as(0)/aho 



FIG. 2: (Color online) Solid lines show the relative energies 
El"^ with / = for the trapped two-boson system interact- 
ing through Vps, Eq. (flT)) with as(_E2°') replaced by as(0), 
as a function of the s-wave scattering length as (0)/aho. The 
dashed and dotted lines show Eq. (|18|) for i < 1 and j = 0, 
and i < 4 and j = 0, respectively. 



TABLE I: Expansion coefScients ^2^1 see Eq. (fTS]) . for 
the weakly-interacting two-boson system with s-wave inter- 
actions. 



* J 



1 

2 

3 

4 



2 1 

3 1 



d 



(ij) 



(_l)"+i23/2/[„!r(-n-l/2)] 

-AK,i/[nlr{-n-l/2)f 

(„l)"+i23/2 {h„,2 + 3hl,) I {n\ V (-n - 1/2)]^ 

(8/3) (/i„,3 + 6/i„,2/i„,i -H 8/1^1) / [n! r (-71 - 1/2)]* 



i(2n + 3/2)4^f 
l(4r)V(2n + 3/2)4^f 



are reported in Table ID The effective range enters first in 
combination with the square of the zero-energy scatter- 
ing length. For the Gaussian potential V^ considered in 
SecHVl the product [as(0)]^re goes to zero as |as(0)| — ?> 
(see solid line in Fig.[T]). Dashed and dotted lines in Fig. [5] 
show Eq. (|18p with j = for i < 1 and i < 4, respec- 
tively. The Taylor expanded expressions for the ground 
state with i <\ and i < 4 agree with the exact eigen en- 
ergy to better than 1% for -0.22 < as(0)/aho < 0.54 and 
—0.67 < as(0)/aho < 0.52, respectively. The accuracy of 
the Taylor expansion deteriorates more quickly for the 
excited states. References [Tfl [13 discuss the structure 
of Eq. (flSj) with i + j < 3 for the ground state as well as 
extensions for N > 2. 

We also expand around the strongly-interacting 
regime. For \ai,olO's[E^2^)\ <^ 1, we expand Eq. (fTT]) 
around the relative energies i?™ 
l/2)huj with n = 



E-f = (2n 



""'' at unitarity, where 
0, 1, • • • . Replacing 



TABLE II: Expansion coefficients di,''^\ see Eq. ^, for the 
strongly-interacting two-boson system with s-wave interac- 
tions. 



^ J 





1 

2 

3 



1 

1 1 

2 1 



2 

1 2 



3 



d 



'(i,j) 



(-l)"+i23/V [2 n! r (-n + 1/2)] 

-4^„,i/[2n!r(-n + l/2)]^ 

(-l)"+i23/2 ^f^^^^^ _^ 3/^2^^ j / [2 n! r {~n + 1/2)]^ 



-i(2n + l/2)<„' 



■(1,0) 




i(4r) +(2n + l/2)4„ 



-|[4:;f^€f + (2n+l/2)4^f 



dh'"r + (2n + l/2)dT„«' 



■(1,0) T(2,0) 



+ 6(2n-h 1/2)4";"^ d. 



+3(2n + l/2)^4-^;^ 



2 7(3,0) 



-i(2n + l/2) 



<2,0) 
•'2,n 



(4:f) +3(2n + i/2)4:;°'4 

+ (2n + l/2)^43f 



l/as(i?2^') by the right hand side of Eq. ([5]), we find 



prcl punit 



j+j<3 



7(M)^«.(0)V7^e V, 



aho / \aho 



^IClI 



7(1,0) 



V 






(20) 



Similarly to the weakly-interacting case, we find it conve- 
nient to express the expansion coefficients ^2 „ in terms 
of the function 



hn.p — Hn^p + (—1) -ff-(i-l/2,j; 



(21) 



Tablellllshows the expansion coefficients dj „ with i+j < 



3. Equation ((20)) contains contributions that are directly 
proportional to the inverse scattering length 1/0^(0), the 
effective range r^ and the volume term. Which of these 
terms dominates depends on the interaction potential. 
Importantly, while the effective range diverges for the 
Gaussian model potential in the |as(0)| — > limit, it re- 
mains finite in the |as(0)|~^ — ^ limit [see dashed line in 
Fig. [IJb)]. The next order terms in Eq. (PH)) are propor- 
tional to [as{0)]-'riV'' with i + j + 3k = A. 

Next we discuss the occupation numbers of the two- 
boson system. The determination of the one-body den- 
sity matrix for the two-boson system requires that the 




Ef/(Iico) 



FIG. 3: (Color online) Occupation numbers per particle 
rivoa/N of the two-boson system interacting through Vps with 
as{E^^) replaced by as{Q) as a function of the relative two- 
body energy El"^ ior Q ^ L ^ M = 0. (a) Solid, dashed 
and dash-dotted lines show nvoo/2 for i/ = 0, 1 and 2, respec- 
tively, on a linear scale. The dash-dot-dotted line shows the 
leading-order depletion, i.e., the first two terms on the right 
hand side of Eq. (f^ . near £2°' = 3&j/2. (b) Lines show 
nvoo/2 for v = Q — &, from top to bottom at £^2'^' — —5huj, 
on a log scale. Arrows mark the local maxima of n^oo/2. 
Near the non-interacting energies, {E2' — E2%)/{huj) is to a 
good approximation directly proportional to the s-wave scat- 
tering length as(0)/aho [see Eq. (|18l) ]. implying that the figure 
can be read as a "scaled occupation number versus scattering 
length" plot. 



wave function tp, written above in terms of the center 
of mass and relative coordinates R12 and ri2 , be trans- 
formed to the single particle coordinates ri and r2. In 
the following, we discuss the occupation numbers asso- 
ciated with the one-body density matrix as a function 
of the relative two-body energy i^j*^', assuming that the 
center of mass wave function is in the ground state, i.e., 
we set Q = L = M = 0. For two-boson systems in one 
dimension, the one-body density matrix was calculated 
as a function of temperature in Ref. |31| . Here, we con- 
sider the zero temperature limit and restrict ourselves 
to relative states with I ~ 0. The formalism developed, 
however, can be straighforwardly applied to states with 
finite I, Q, L and M . As detailed in Appendix |B1 the 
one-body density matrix can be evaluated efficiently and 
with high accuracy by expanding it in terms of single 
particle harmonic oscillator functions. 

We first consider the zero-range pseudopotential Vps 
with as{E^2^) replaced by as(0). Figure [3] shows the 
scaled occupation numbers n^oo/2 for the trapped two- 
boson system interacting through V^s as a function of the 
relative two-body energy i?2°'- In the non- interacting 



limit, the ground state with energy £2^^ = 3Sw/2 is 
characterized by a single natural orbital with projec- 
tion A = 0, i.e., the non-interacting two-boson system 
in the ground state has a condensate fraction N^/N of 1. 
As the interactions are turned on, i.e., as 0^(0) takes 
on small positive or negative values corresponding to 
E2°^ > 3tiuj/2 and E2°^ < 3haj/2, respectively, the occu- 
pation number associated with the lowest natural orbital 
depletes. Taylor-expanding the projected one-body den- 
sity matrix around £'2°' = 3fta;/2 (see Appendix |B|) . we 
find that the condensate fraction of the ground state of 
the two-boson system depletes quadratically with as(0). 



n = 0]. 



No/2 



1 - 0.420004[a^(0)/aho]^ 
-0. 373241 [a^(0)/aho]^ 
+ 0.406786[a^(0)/aho]'' 
+ O(K(0)/aho]'). 



(22) 



The dash-dot-dotted line in Fig. [2Ia) shows the first two 
terms of Eq. (|22l) . i.e., the leading-order dependence of 
the depletion on 0^(0), in the weakly- interacting regime. 
The higher-order corrections proportional to [0^(0)]% i = 
3 and 4, are analyzed in Sec. |lVl 

Figure [3] reveals oscillatory behavior of the scaled oc- 
cupation numbers n^oof^ for £2'^^ > ihjj/2. As the rel- 
ative energy E^2^ {E'^'^ > 3huj/2) increases, the occu- 
pation numbers go through "near deaths and revivals" , 
with many of the occupation numbers crossing. When 
a higher-lying non-interacting state is reached, one more 
natural orbital with A = becomes macroscopically oc- 
cupied. Similar structure is seen for A > (not shown). 
For E2°^ — 7huj/2, e.g., five natural orbitals are occu- 
pied. Two of these natural orbitals have A = with 
ni/oo/2 = 1/4 (solid and dashed lines in Fig. [3|), while 
three natural orbitals (from the 2A -|- 1 degeneracy) have 
A = 1 with n,yoo/2 = 1/6. The largest occupation num- 
ber per particle nooo/2 takes on local maxima at the non- 
interacting energies E2°^ = llhjj/2 and 19fta;/2 as well 
as at £'2°' ~ 3.68?icj and 7 .QAhuj (see arrows in Fig. [3]). 

For £2°' < 3?iaj/2, the scaled occupation number 
^000/2 (see solid line in Fig. [3]) decreases monotonically 
with decreasing energy while many other natural or- 
bitals become occupied, including natural orbitals with 
A > 0. In the limit of a deeply bound two-body state, 
the relative two-body wave function becomes infinitely 
sharply peaked, implying that infinitely many single- 
particle states are required to describe the deeply-bound 
two-boson system. 

If the energy-dependence of the s-wave scattering 
length is accounted for, the condensate fraction of 
the two-boson system in the ground state near the 
non-interacting regime depends not only on 0^(0), see 
Eq. (|22|) . but also on r^- We find that the leading-order 
eff'ective range contribution to the condensate fraction is 
-(3/2) X 0.420004re[a^(0)] Vajio- The factor of 3/2 arises 
since the scattering length as(0) has to be replaced, fol- 
lowing Eq. ([6]), by [as(0)]^refc^/2, where the relevant en- 
ergy scale for evaluating k^ is ifujj/2 [see Eq. ^S]) for 



IV. 



WEAKLY-INTERACTING TRAPPED 
iV-BOSON GAS 



Subsection IIV Al determines the condensate fraction of 
the lowest gas-like state of the A^-boson system perturba- 
tively in the weakly-interacting regime, and compares the 
perturbative predictions with our results for finite-range 
potentials. Subsection IIVBI parametrizes and quantifies 
the non-universal corrections revealed through the com- 
parison. 



A. Perturbative treatment 

We employ the formalism of second quantization and 
rewrite the A'^-boson Hamiltonian as 

iJ = ^ E'aat fla + ■;^ X! -f^^abcdalabOdac, (23) 



abed 



where 

ATabcd = 



(24) 



^$:(ri)$£(r2)14b(ri - r2)$c(ri)$d(r2)d'rid3r2. 

Here, the $a(r) denote the single particle harmonic os- 
cillator wave functions with eigen energy E^.] in spherical 
coordinates, we have Eg^ = {2na + la + 3/2)huj. The oper- 
ators Oa and a^ obey bosonic commutation relations and 
respectively annihilate and create a boson in the single 
particle state $a- We model the interaction Vtb by the 
sum Vp + V' (see Se c.lnl) , and employ the counterterms 
derived in Refs. [1^ [l3| to cure divergencies. Since the 
<I'a(r) are known, the matrix elements A'abcd for this in- 
teraction model can be evaluated either analytically or 
numerically [16l [l7| . The low-energy Hamiltonian given 
in Eq. (j23p has previously been used to derive pertur- 
bative energy expressions up to order a^^ for the har- 
monically trapped A^-boson system [l3|. In particular, 
Vp and V were treated at the level of third- and first- 
order perturbation theory. Comparison with energies for 
systems with finite-range interactions validated the for- 
malism and showed that the derived perturbative energy 
expressions, which depend on as(0) and rg, provide an ex- 
cellent description in the weakly- interacting regime [17 1. 
To determine the condensate fraction, we construct the 
matrix (aj^aq), where p and q run over all possible single- 



particle state labels. The expectation value (ajjOq) 



(flpaq) 






(25) 



is calculated with respect to the many-body ground state 
wave function i/jq , determined within fc^'^-order pertur- 

(k) 

bation theory. The ground state wave function ■0o is 



expressed as a superposition of the unperturbed many- 



body wave functions ip\ 



(0) 



1^^' 



(fc)x 



E^r 



i^f) 



dk) 



(26) 



where the expansion coefficients b- are determined by 
the matrix elements if abed- The subscript j coUectively 
labels the non-interacting or unperturbed many-body 
states (in particular, labels the ground state). Once 
the matrix {a^ciq) is constructed, we diagonalize it ana- 
lytically (see Appendix R)). Up to order a^^, i.e., treating 
the second-order perturbation theory wave function, we 
find 



a,(0) 



1 2 



Oho 



No/N = 1 - 0.420004(iV - 1) 

+ [-0.373241(iV-l) 
-h0.439464(A^-l)(iV 



as(0) 



aiio 



(27) 



The prefactors are discussed in the context of Table IIIII 
We interpret the terms proportional to (iV — 1) and 
{N — 1){N ~ 2) as being due to two-body and three- 
body scattering processes, respectively. In Eq. (|27p . the 
terms proportional to [as(0)]^ and [as(O)]"^ arise when 
treating the potential Vp, together with the appropri- 
ate counterterm, in first- and second-order perturbation 
theory. Treating Vp in third-order perturbation theory 
(not pursued here), three terms proportional to [as(O)]'* 
that contain the factors (TV - 1), {N - 1){N - 2) and 
{N — 1){N — 2){N — 3), respectively, are expected to 
arise. 

The potential V' does not, in first-order perturba- 
tion theory, give rise to a two-body term proportional to 
re[as(0)]^/a^o- This result agrees with that obtained by 
Taylor-expanding the full two-body density matrix and 
determining its largest eigen value (see last paragraph of 
Sec. IIII|) . No three-body term arises at order r'e[as(0)]'^. 
Since the leading-order effective range dependence is of 



order a^^ , it is not included in Eq. (j27p . 

To assess the applicability of our perturbatively de- 
rived result, Eq. ([27]), we calculate the condensate frac- 
tion for small A^-boson systems interacting through 
the Gaussian model potential Vg, Eq. ([TO)) . with small 
|as(0)|/aho- For A = 2, we solve the relative Schrodinger 
equation using standard B-spline techniques. For A = 3 
and 4, we use the stochastic variational approach [32|,|33|, 
which expands the relative eigen functions in terms of a 
set of fully symmetrized basis functions whose widths 
are chosen semi-stochastically. The widths are optimized 
by minimizing the ground state energy. The optimized 
ground state wave function is then used to construct the 
projected one-body density matrix poo(^ij^i) on a grid 
in the r'l and ri coordinates. The projected one-body 
density matrix is diagonalized numerically to find the 
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FIG. 4: (Color online) Condensate fraction Nq/N of the 
weakly-interacting A^-boson system as a function of Os(0)/aho. 
(a) Squares and circles show the condensate fraction for N = 2 
and 3 bosons interacting through the shape-dependent Gaus- 
sian model potential Vg with ro — O.Olaho. Solid and dashed 
lines show Eq. (|27p for N = 2 and 3, respectively, (b) Squares 
and circles show the quantity No/N - (Nq/N)^^^ ior N = 2 
and 3 using the same data as in panel (a). Solid and dashed 
lines show the [as(0)/aho]^ term of Eq. ^7} for iV = 2 and 3. 



natural orbitals and their occupation numbers. The re- 
sulting condensate fraction for A^ = 3 has an estimated 
numerical error of order 10^^ or smaller for the parame- 
ter combinations considered. The numerical error is due 
to the facts that (i) we use a finite basis set, (ii) we con- 
struct pQo{r[, ri) on a grid with finite grid spacings, and 
(iii) our grid terminates at finite r'l and ri values. For 
A = 3, we use around 120 basis functions and 625 lin- 
early spaced grid points in r[ and ri between ro/2 and 
5.5aho. 

Squares and circles in Fig. Ufa) show the condensate 
fraction Nq/N for the N — 2 and A = 3 systems inter- 
acting through Vg with tq — O.Olaho as a function of the 
zero-energy s-wave scattering length as(0). For compari- 
son, solid and dashed lines show our perturbative results 
up to order a'^^ . It can be seen that the numerically de- 
termined condensate depletions for the finite-range inter- 
action potential change approximatically quadratically 
with the scattering length. To investigate the correc- 



tion proportional to [as(0)/aho]^, squares and circles in 
Fig. |li;b) show the quantity Nq/N - {No/N)'^'^\ where 
{No/N)'^^'> = 1 - (iV- l)0.420004[a,(0)/aho]^ for TV = 2 
and 3 [the data for Nq/N are the same as those shown in 
Fig. m^a)]. For comparison, the solid and dashed lines 
show the perturbatively predicted [as(0)/aho]'^ depen- 
dence for TV = 2 and 3. It can be seen that the per- 
turbative prediction does not provide a good description 
of the sub-leading dependence of the depletion. Similar 
behavior is observed for iV = 4 (not shown) . We find that 
neither the inclusion of the rg[as{0)]^ nor of the [0^(0)]^ 
terms can explain the discrepancy. As shown in the next 
subsection, the discrepancy displayed in Fig. \Mjo) is due 
to non-universal contributions. 



with r'i2 = r'l - r2 and R12 = (r'^ -I- r2)/2; ri2 and R12 
are defined in Sec. |IT1 In writing Eq. (P^. the term pro- 
portional to \Stp\'^ has been neglected. 

To determine the condensate fraction, we expand 
Splr'i, ri) in terms of non-interacting single particle har- 
monic oscillator functions; this approach is analogous 
to that discussed in detail in Appendix IB] for p{r[,ri). 
A fairly straightforward analysis shows that the main 
correction to the condensate fraction arises from the 
((0,0,0), (0,0,0)) element of Sp{r[,ri). It follows that 
the largest occupation number nQQQ for the finite-range 
potential can be written as 



'000 



'^ooo 



+ Scoo, 



(31) 



where 



B. Non- universal contributions 

To better connect the results obtained by applying per- 
turbation theory to the low-energy Hamiltonian, Eq. (j23p 
with I4b = Vp + V' , with the results for the Hamiltonian 
with finite-range interactions, we first consider the two- 
body system. Noting that the perturbative results for the 
condensate fraction agree with the results obtained by 
Taylor-expanding the exact one-body density matrix for 
the two-body system with regularized zero-range inter- 
action [i.e., noting that Eq. (P7)) agrees with Eq. (|22p for 
N = 2 and the orders considered], we compare the one- 
body density matrix p(r'i, ri) for the two-body system in- 
teracting through V^s with the one-body density matrix 
/9fr(r']^, ri) for the two-body system interacting through a 
finite-range potential. 

Since the regularized zero-range potential reproduces 
the relative two-body energy of systems with finite-range 
interactions with high accuracy [l8|, |l9| , we assume that 
the relative two-body energies £2^^ agree for the two 
interaction models. We denote the normalized relative 
wave function of the energetically lowest-lying gas-like 
state for the finite-range potential by ■(/'"^''^'(ri2) and 
that for the regularized zero-range model potential by 
■0^00 (''12) [see Eq. (|B2I) ]. It is instructive to write ijj'^'^^'^'' 
as 



V.-i.f(ri2)=V^JSo(ri2)+<5VXri2). 



(28) 



Inserting Eq. (^5)) into pfr{r[,ri) and assuming the ab- 
sence of center of mass excitations, we find 

Pft.(rl,ri) «p(rl,ri)+2(5p(r;,ri), (29) 

where p{r[, ri) is given in Eqs. (jBip and (IB16|) . and 

Mr'i,ri)- (30) 

[^^{ro(R'i2)^^(r'i2)]*«o(Ri2Xoo(ri2)rf=^r2 

+ /[^SSMR'i2)V'roo(r'i2)]*«o(Ri2)^^(ri2)d^r2 



fcoo = 5l2"''(CrA + C,A*)- 



(32) 



In Eq. pil) . nooo/2 — Nq/2 denotes the condensate frac- 
tion of the two-body system interacting through the reg- 
ularized zero-range potential [see Eq. (|22p ]. The coeffi- 
cients Ci are defined in Eq. (|B8I) and the Di denote the 
overlaps between the non-interacting harmonic oscillator 
functions and S^^ 



D, 



= / [^IS!,^"'(ri2)]*<SV(ri2)d^ri2. 



(33) 



Realizing that the z = terms in Eq. (|32p dominate and 
using the leading-order behavior of Co, i.e., Co ~ 1 [see 
Eqs. (|BT9l) and (fB20| ]. we find 



,5coo « 2Re{Do). 



(34) 



For the finite-range potentials considered in this paper, 
we find that Eqs. ([32]) and ^ deviate by less than 0.2 %. 
In the following, we refer to Dq as the "non-universal 
two-body parameter" . 

Figure [5] illustrates the behavior of the integrand that 
determines the non-universal two-body parameter Dq for 
two different two-body energies, i.e., for £^2°' — l-498?ia; 
[negative as(0)] and i?2°' = 1.502huj [positive 0^(0)], for 
the Gaussian model potential with vq — O.Oloho- In 
particular. Fig. [S{a) shows the product (iPq'qq^^)* Sip and 
Fig.[5jb) the ratio V''°^'^'/'0noo- Figure [5] reflects the pres- 
ence of the two characteristic length scales of the prob- 
lem. The behavior in the small ri2 region, shown in the 
insets of Figs.ISJa) andlHJb), is governed by the details of 
the two-body interaction potential. Near ri2 ~ 5ro, the 
behavior of (V'ooo"')*'^^ ^^'^ ^''^''''^'^/V'qoo changes notably. 
For ri2 > 5ro, the ratio '^'^''''^'^/V'goo approaches a con- 
stant that is slightly larger (smaller) than 1 for negative 
(positive) as(0). The small deviations of the ratios from 
one are a consequence of the fact that the wave func- 
tions of the trapped system are normalized to one. Since 
the wave functions for the finite range interaction po- 
tential deviate from those for the zero-range interaction 
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FIG. 5: (Color online) Analysis of the integrand that deter- 
mines the non-universal parameter Do [see Eq. (|33p ] for the 
finite-range model potential Vg with ro = O.Olaho for two dif- 
ferent relative two-body energies. The solid and dotted lines 
correspond to i?^' = 1.498&J and Eq"^ = 1.502/ia;, respec- 
tively. Panel (a) shows the integrand ['(/)QQQ'"(ri2)]*(5'(/'(ri2) in 
the large ri2 region, i.e., for ri2 G [O.OSaho, 3aho] (log-linear 
scale), while the inset shows the absolute value of the in- 
tegrand |[i/'ooo"'(''i2)]*5i/;(ri2)| in the small ri2 region (log- 
log scale). Panel (b) and the inset of (b) show the ratio 
^roi,fr^^roi^ in the large and small ri2 regions, respectively 
(log- linear scale). 



potential in the small ri2 region, the ratio '/''^"''^'^/V'qoo 
needs — in general — to differ from one in the large ri2 
region. The "divergence" of the ratio ■i/''^''''^'^/'0qoo near 
ri2 = 0.002aho for 0^(0) > [see dotted line in the inset 
of Fig.ISjb)] reflects the fact that the zero-range potential 
supports a deeply-bound negative energy state, which in- 
troduces a node at small ri2 in the wave function that 
describes the energetically lowest-lying gas-like state. A 
corresponding bound-state is not supported by the purely 
repulsive Gaussian model potential, leading to an infinite 



at the node of V'Joo ■ 



ratio ^'^^•^'/^gfo 

The non-universal two-body parameter Dq is the trap 
analog of the two-body scatteringquantity uq introduced 
by Tan [see Eq. (114a) of Ref. Q]- It is important to 
note, however, that Dq depends on the wave function 
difference for all ri2 and on the non-interacting harmonic 
oscillator ground state wave function while the uq defined 
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FIG. 6: (Color online) Squares and circles show the residual 
No/N - (N/No)'-^^ for iV = 2 and 3 as a function of as{0) for 
the Gaussian model interaction Ig with ro = O.Olaho- Solid 
and dashed lines show the quantity (A'' — l)5coo for N = 2 
and 3, respectively [3J|. 



by Tan depends only on the wave function difference in 
the small ri2 region, i.e., out to a few times tq. Indeed, 
we find that the contribution to Dq that accumulates 
in the inner region (ri2 ^ lOro) can be of comparable 
magnitude to the contribution that accumulates in the 
outer region (ri2 > lOrg). 

Since the dependence of the condensate fraction on the 
quantity Dq arises at the two-body level, the Dq term 
needs to be multiplied by TV — 1 for systems with N > 2. 
Figure [6] compares the condensate fractions for the N = 2 
and 3 systems interacting through the finite-range Gaus- 
sian model potential Vg with the predicted behavior for 
the condensate fraction. Specifically, squares and circles 
show the difference Nq/N - (Nq/NY^^ for iV = 2 and 3 
between the numerically determined condensate fraction 
N/Nq and the perturbative result {Nq/NY^\ which in- 
cludes all terms on the right hand side of Eq. ([27)) up to 
order a^^. According to our discussion above, we expect 
that the residuals are well approximated by {N — 1)Scqq 
(shown by solid and dashed lines in Fig. [5]). Indeed, Fig.[S] 
shows that the residuals for A^ = 2 and 3 are well de- 
scribed by the non-universal corrections. We find similar 
results for A^ = 4 (not shown). Our calculations demon- 
strate that two Hamiltonians that are characterized by 
the same energy give rise to condensate fractions that dif- 
fer. Related findings have previously been discussed in 
Refs. p, 35, 36]. The leading order difference between the 
condensate fractions for the harmonically trapped few- 
boson systems described by the two Hamiltonians can be 
parameterized by the non-universal two-body parameter 
Dq, a parameter not needed to match the energies of the 
two Hamiltonians. 



THREE TRAPPED BOSONS AT UNITARITY 
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The previous section discussed the condensate fraction 
of weakly-interacting trapped iV-boson systems, which 
can be expressed in terms of as(0), r,, and a non- universal 
two-body correction parametrized through Dq. It is well 
known that the properties of the three-boson system 
not only depend on two-body parameters, but also on a 
three-body parameter |37H39| . In the weakly- interacting 
regime, however, the dependence on the three-body 
parameter appears at higher order than considered in 
Sec. IIVI In the strongly-interacting regime, in contrast, 
the dependence on the three-body parameter is generally 
quite pronounced. 

At unitarity, i.e., for diverging as(0), the trapped 
three-boson system with zero-range s-wave interactions 
supports two distinct classes of eigen states: (i) universal 
states whose properties are fully governed by the two- 
body scattering parameters, and (ii) non-universal states 
whose properties depend, in addition to the two-body 
scattering parameters, on a three-body parameter. In the 
following, we determine the occupation numbers for the 
non-universal three-boson states in a trap at unitarity as 
a function of the three-body parameter. The momentum 
distribution of Efimov trimers in free-space was discussed 
inRef. ^. 

The three-boson wave function V'lri, r2, '^3) with rela- 
tive orbital angular momentum / = for zero-range in- 
teractions with diverging s-wave scattering length as(0), 
and vanishing rg and V , under external isotropic har- 
monic confinement can be written as [4Cl l4l| 

V'(ri,r2,r3) =5 [i?-5/2^(i?)^(a)V'^'^M(Ri23)] • (35) 

Here, R denotes the hyperradius and a the hyperan- 
gle, R^ = rl^/2 + 2rl2^j2, and tana = %/3ri2/(2ri2,3) 
with ri2 = |ri - r2| and ri2,3 = |(ri + r2)/2 - raj. 
In Eq. ([55]) . ■!/'ni,M(^i23) denotes the harmonic oscil- 
lator wave function in the center of mass coordinate 
R123, R-123 = (ri + r2 + r3)/3; as in Sec. [nil we as- 
sume that the center of mass wave function is in the 
ground state, i.e., we set Q = L = M = 0. The op- 
erator S ensures that the three-boson wave function is 
symmetric under the exchange of any of the three bo- 
son pairs, S = 1 + P12+P23+ P31 + P12P23 + P12P31, 
where Pjk is the operator that exchanges particles j 
and k. The hyperangular wave function (p{a) takes 
the form ip{a) = sin [(a — 7r/2)so]/sin(2a), where sg 
equals 1.00624i |37H39l |. The fact that the separation 
constant so, which arises when solving the hyperangular 
Schrodinger equation, is imaginary is unique to the I = 
channel and tightly linked to the fact that a three-body 
parameter is needed. 

The hyperradial wave function F{R) can be conve- 
niently expressed in terms of the Whittaker function 
W m, i-e., F{R) = i?-i/2W^B-/2,.o/2(^V«L)- The 
relative three-body energy E'™' is related to the three- 
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FIG. 7: (Color online) Occupation numbers per particle 
n-voo/N [i' = (solid line) and u = 1 (dashed line)] for the 
non-universal I = state of the three-boson system at unitar- 
ity as a function of the relative three-body energy £^3''' in units 
of the oscillator energy -Eho, i?ho = huj; the three-body energy 
can be converted to the three-body parameter via Eq. 



body or Efimov phase 6 through J^C 

f- 
= arg' ^^ 




(36) 



The physical meaning of 9 becomes clear when look- 
ing at the small R/aho behavior of F{R), F{R) — > 
vi?sin(Im(so) \n{R/aho)+S). This expression shows that 
the three-body phase determines what happens when 
three particles come close together. The smaU R/uho 
behavior can be thought of as being imposed by a short- 
range three-body force or a boundary condition of the 
hyperradial wave function in the i?/aho -^ limit 39] . 

To determine the occupation numbers of the non- 
universal three-boson states as a function of E^°^, we sam- 
ple the density [-(/"(ri, rg, rs)^ using Metropolis sampling. 
As discussed in Ref. [28] , this approach introduces a sta- 
tistical error that can be reduced by performing longer 
random walks. Throughout our random walk, we sample 
the projected one-body density matrix poo(^ii^i)- Diag- 
onalizing pooiTi, ri) at the end of a run yields the occu- 
pation numbers. 

Figure [7] shows the two largest occupation numbers 
per particle ni,oo/3 as a function of E^°^. It can be seen 
that the occupation numbers of the non-universal state 
depend quite strongly on the relative three-body energy 
or, equivalently, the three-body phase 9. The maximum 
of the lowest occupation number per particle nooo/3 is 
0.82 and occurs at £^3''' = 3huj/2; the occupation num- 
ber per particle nioo/3 is minimal at this energy. In- 
terestingly, the occupation numbers show oscillations (or 
"shoulders" ) similar to those discussed in the context of 
Fig. [3] for the two-body system. In Figure [3l we change 
the relative two-body energy, which is related to the s- 
wave scattering length through Eq. (IT71) . In Fig. [71 we 
change the relative three-body energy, which is related 
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to the three-body phase through Eq. ((36|) . For both the 
two- and three-body systems, Eqs. ([T7]) and ([55]) can be 
related to the short-range boundary condition of the re- 
spective radial or hyperradial part of the relative wave 
function. 

For comparison, we also calculated the largest occupa- 
tion number per particle of the projected one-body den- 
sity matrix popfr ^ , ri ) for selected universal three-boson 
states (see Ref. [41j for the relevant wave functions). 
The largest occupation number per particle of the en- 
ergetically lowest-lying universal three-boson state with 
Z = and ^J°^ = 4.465fia; at unitarity is nooo/3 = 0.295. 
The largest occupation numbers per particle of the en- 
ergetically lowest lying states with / = 1 and / = 2 
at unitarity are nooo/3 = 0.199 and 0.421, respectively. 
For these states, the energies are E^'^^ = 2.864fiLu and 
El°^ = 2.823fiaj, respectively. For the universal states 
considered, the largest occupation number riooo is no- 
tably smaller than N. 



VI. CONCLUSIONS 

We have determined and interpreted the occupation 
numbers of few-boson systems under isotropic harmonic 
confinement. In the weakly-interacting regime, our anal- 
ysis is based on a low-energy Hamiltonian — characterized 
by the s-wave scattering length 0^(0) and the effective 
range r^ — that has previously been proven to correctly 
describe the energy of few-boson systems up to order 
«ho UM- The present paper shows that this low-energy 
Hamiltonian correctly describes the leading order deple- 
tion of harmonically trapped few-boson systems but that 
it does not fully capture the corrections to the leading 
order depletion. 

Our final expression for the condensate fraction reads 



a.(0) 



flho 



Nq/N = 1 - 0.420004(iV - 1) 

+ [-0.373241(7V-1) 

+ 0.439464(iV-l)(iV-2)] 



.(0) 



flho 



+ [0.406786(iV-l) + 7_f^(A^-l)(7V-2) 



+ ji'^'^iN - 1){N - 2){N - i)] 
2{N-l)Rc{Do) 



as{0) 



Oho 



(3/2) X 0.420004(iV- 1 



reiasiO)] 



^ho 



(37) 



where the non-universal two-body parameter Dq is de- 
fined in Eq. ((33|) . The coefficients 73 and 74 arise when 
treating Vp in third-order perturbation theory; the deter- 
mination of their numerical values is beyond the scope of 
this paper. We have confirmed the expression for the 



condensate fraction, Eq. (|37| . through comparison with 
numerical results for a few-body Hamiltonian with finite- 
range two-body potentials. Our work demonstrates that 
the occupation numbers are not fully determined by the 
parameters of the "usual" effective range expansion, but 
rather depend on an additional property of the two-body 
wave function (i.e., non-universal physics). A similar re- 
sult is expected to hold for the momentum distribution. 
Our findings are not only of importance for cold atomic 
Bose gases but also for nuclear systems, for which the 
use of low-energy Hamiltonians has become increasingly 
more popular during the past decade or so |43j . 

We have also considered the strongly-interacting 
regime. Our results show that the occupation numbers 
for non-universal states of the three-boson system un- 
der isotropic harmonic confinement depend strongly on 
the three-body parameter. This finding suggests that 
the occupation numbers and momentum distribution of 
strongly-interacting Bose gases at unitarity may depend 
on three-body physics. In view of recent experimental 
work pil - liq , it would be interesting to extend the treat- 
ments of Refs. [43, |4^ , which predict — accounting only 
for two-body physics — that three-dimensional Bose gases 
at unitarity fermionize. In particular, it would be in- 
teresting to determine how, if at all, this fermionization 
picture changes if three-body physics is accounted for. 
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Appendix A: Perturbative treatment and 
diagonalization of one-body density matrix 

This appendix provides details regarding the perturba- 
tive treatment of the condensate fraction of the iV-boson 
system and the diagonalization of the associated matrix. 

We start with the Hamiltonian H given by Eq. (|23p . 
We first neglect the effective range dependent potential 
V , i.e., we consider only the bare Fermi pseudopotential 
Vp, Eq. ([5|, and the counterterm W used to cure diver- 
gencies [la, II4I • The matrix elements i^abcd can then be 
written as 



K, 



abed 



F, 



abed 



aho 



flho 



fiw, (Al) 



where 



abed 



4™^ I <l>:(ri)$* (ri)<i>e(ri)<i>d(ri)d3ri (A2) 



and 



t.= \/-(l-ln2) 



2^^ 



(A3) 
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The coefficient /^j was calculated in Ref. [l3l and is 
listed in Table IIIII It diverges and cures the divergencies 
that arise when treating Vp in second-order perturbation 
theory. 

is evaluated by substituting 

In order to get the matrix ele- 

we need to employ second order 



(dta„ 



The matrix 
Eq. daSl) into Eq! ^'. 
mcnts up to order a^^ . 

perturbation theory. The expansion coefficients b^' of 

(2) 

the non-normalized second-order wave function ^q read 



(2) 



^(2) 



(A4) 



and 



for j not equal to the ground state labeled by 0. In the 
denominators appearing in Eq. (|A5|) . the E- denote the 
unperturbed eigen energies corresponding to the j's un- 
perturbed eigen state. The numerators are conveniently 
expressed in terms of the matrix elements i^abcd • 

The indices p and q of (a^aq) run over all possible sin- 
gle particle state labels. We employ spherical coordinates 
and write p — {n'^,l[,m'^) and q = (ni,/i,mi). We find 
that the matrix is block diagonal, i.e., (a]^, ;, ^_^^, an-^i-^m-^) — 
for l[ ^ li or m[ ^ rrii. In the following, we consider 
the submatrix with I'l = /i = m'l = nii = 0. We denote 
the matrix elements by c„' m and write 



b^^^ 



do) 



{4'\'"\Vf + W\^'^ 



?{o) 



(0)\ 



E^'-Eo 



(0) 



■E 

-,(0) 



(tA "Vf + w^i4"0(4 Vf + w^i^!," 



(E: 



(0) 



(o)\/,/,(o)| 



4")) 






',(0)\ 



(A5) 



C<n.=j:-^ny+0{x% 



fc=0 



(A6) 



where 2: = as{0)/aho- Considering symmetry and keeping 
terms up to order a;^, we find 



"■n'OO^-mOo) 



/ 1 

(1) 



(2) 
'^00 ' 



(3) 



10 10 10 



(1) 



(2) 2 



-10 •' 



„(3)„3 



(2) 
-10 ■" 



„(3). 



-10 ■ 

(2) ,_^ (3)^3 



^(2)„2 , „(3) 3 



„(1), 



J2)^2 



(3)^3 \ 
JnX 



-AO 



(2) ,_^ (3)^3 



^(2) 2 



-Al- 



■ C A aX 



-AA-' 



(A7) 



The upper left element is 1, with small corrections pro- 
portional to x^ and x'^. The leading-order contribution 
of the other elements in the first row and first column is 
proportional to x. The leading-order contribution of the 
rest of the matrix elements is proportional to x^ . 
We diagonalize the matrix by solving 



det(M) = det((a,V 



oo^moo) 



SI) 







(A8) 



through application of the Leibniz formula for determi- 
nants f49!|. In Eq. (|Mt . I denotes the (A + l) x {A+ 1) 
identity matrix and S the eigen value we are seeking. The 
product of the diagonal elements can be written as 



by Mik and Mki- For fc = 2, for example, we have 



M12M2 
M11M22 



A+l 



- n ^^" 



{^x+^x^+^x^yiK^x^+c 



(2)^2^ J3)^3_2 



i=2 



(c«)^x^ + 2c«cfo)x3l i-E)^-' + O (x^) . (AlO) 



Summing over all contributions with fc > 2, we find 



A+l 


=(-s)-4+i + 


A 


i=l 




L i=o 




A 







^^^-.4^3 



"\A 



(-S) 



x"|(-S)^-i+0(x4) 



(A9) 

The other terms involve the product of the diagonal ele- 
ments with the first and fc*'' diagonal elements replaced 






E \i^ 



))^x^+2c«cf)x3 



"-jo >-io 



i=i 



(-S)^-i-HO(x4). 

(All) 

Combining Eqs. (jA9[) and (IA11[) yields the eigen value 



13 



equation up to order x^ 

A 



(-5) 



77\A+1 



3=0 



(....,(3)^3 



(-S)^ 



E 



J2) 2 , J3)^3 



„(1)^2 2 



(-S)^-i = 0. (A12) 






Equation (|A12p can be reduced to a quadratic equation in 
S. Taking A to infinity, the largest eigen value coincides 
with the condensate fraction, 



No/N = 1 + 



.(2 



l)^2 



00 "T A^^So ) 

.(3) , V^9J1)J2) 



(A13) 



x3 + O (x^) 



The coefficients cLn are determined by Eqs. ([25|). ([26| . 
(IA4[) . and (jA5[) . and can be expressed in terms of infinite 
sums involving the matrix elements -Fabcd (see Table Hill) . 
Evaluating the coefficients cLn, Eq. (JAISP becomes 



(2), 



iVo/A^-l-7r(^-l) 



„(3) 



a,(0) 



flho 



l^iN-l) 



v(3) 



+7r(^-i)(^-2) 



a,(0) 



n3 



ail, 



(A14) 



where 72 = 



-27^^^ - 



■^ifl 



2^^il and 7^ ^ 



-47^ 



The superscript and the first sub- 



script of the coefficient 7^ denote respectively the orders 
of as(0)/aho and the multi-body scattering process that 
7j ■ is associated with. The second subscript simply la- 
bels the various sums (see Table IIIII) . To evaluate 73 4 , 
we use the expression 



73,4 — 



3/2 



00 00 

EE 
j=i fc=i 



TT^ 1 , 

^+ln2--ln^2 



ln(8 - 4\/3) - 

2l/2-2j-2/cp(^- ^ ^. ^ 3/2) 



pkn'^jlkl 



(A15) 



(fc) 



If we insert the numerical values of the coefficients 7^ 
from Table Hill we obtain Eq. (|27p of the main text. 

To understand how the effective range contributes to 
the depletion of the condensate fraction, we treat the 



TABLE III: Expressions for and numerical values of the coef- 
ficients % that enter into Eq. HA14[I . The representation of 

the % ! in terms of infinite sums, derived within the pertur- 
bative frameworlc, are listed in column 2. For completeness, 
we also list the coefiicient fi\ , which enters into the countert- 
erm W needed to cure the divergencies arising from Vf. Atab 
denotes a dimensionless energy; in spherical coordinates, we 
have Aeab = 2na -I- ^a -I- 2n6 -I- li,. The sums are over all vec- 
tor indices with the restrictions a 7^ 0, b 7^ and c 7^ 
(e.g., the sum that determines 72 is X^ = Sa^^o b=.^o' where 
a = corresponds to Ua = la = 0). The numerical values 
for the coefficients are given in column 3: 72 j 72 1, 72 21 72 3 
and 7;^^] are obtained by evaluating Eqs. (IB29|) . (IB30|) . p31|) . 

(|B32[) . and (|A15|) while 73 3 is obtained by evaluating the in- 
finite sum numerically (the numerical uncertainty is reported 
in round brackets). In terms of the a coefficients defined in 



Ref. 


m 


we 


have7f = (4^^-2af 


')M^ 






infinite sum 


numerical value 


P^'' 




E 


FoOab-FbaOO , V ^OOOafaOOO, 
Ae„b ' "^ Ae,o 


diverges 


4'^ 




0.420004291120 


<l 


Foooo E """"g^fi-"" 


0.073250101788 


(3) 

72,2 


y^ f'oOOaf'aOObfbOOO, 


0.005269765990 


(3) 
72.3 


(l-ln2)y^|7f 


0.102830963978 


73,1 


72.1 


0.073250101788 


73,2 


72.2 


0.005269765990 


73,3 


Y^ ^OOabJ^bOOc^caOO 


0.067074(1) 


'3,4 


y^ foOOafaObcJ'cbOO 


0.054238116273 



potential V in first-order perturbation theory. We find 
that the V' does not give rise to a term proportional to 

re[a,(0)]V4o- 



Appendix B: Determination of one-body density 
matrix for N — 2 

This appendix summarizes the evaluation of the one- 
body density matrix for the two-boson system with regu- 
larized (5-function interaction in a spherically symmetric 
harmonic trap. We start with Eq. (fTTll and write the 
two-body wave function as a product of the center-of- 
mass wave function V'qlm (-^^12) ^.nd the relative wave 
function ^p"ll^{rl2)■ In the following, we assume that the 
two-body wave function is normalized and restrict our- 
selves to states with Q — L — M = l = m = 0, yielding 



p(r;,ri) ^2l[roZiK2Kooir^'i2)r 



V'SS^„(Ri2)V'roo(ri2)d'r2. 



(Bl) 



To evaluate Eq. (jBip . wc follow a three-step process: 
(i) We expand the relative wave function in terms of a 
complete set of non-interacting harmonic oscillator wave 
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functions in the relative coordinates, (ii) We expand the 
non-interacting relative and center of mass wave func- 
tions in terms of non-interacting single particle harmonic 
oscillator wave functions, (iii) We integrate over r2. 
Step (i): The relative wave function reads [29i] 






Aw 



.„ 3 1 

y' 2' 2 



ri2 

Oho 



1 f 2112-\'^ 



£12. V 



(B2) 



where U is the confluent hypergeometric function and the 
normalization constant A^™' is given by 



^.ei_,/ 2^^ r(-l - 2g) a->-V^ 23/2 
" ]] l/q + TT cot(7rg) - V(-9 - 1/2) + ^{q) ' 

Here, ip is the digamma function and the non-integer 
quantum number q is determined by the s-wave scatter- 
ing length via Eqs. (TTBl) and ([T7| . In the non- interacting 
limit, we have 



/rel.ni/ \ 
V-iOO (l"l2) 



Nr 



-(1/2) 1 



47r 



ri2 



flho 



e 4 V ai„ 



(B4) 



relative and center of mass wave functions in terms of 
single particle states, 

Em™ ^QiM(Ri2)C™">i2)(i, M, I, m|A, H) = 
J2mhrm E„2i2™2 (("1' 'i' "2, ^2; A|Q, L, i, /; A» x 

(?1, TOi, /2, m2| A, n)$„jii™i (ri)$n2i2m2 (i"2), (BIO) 

where ((•••)) denotes a Talmi-Moshinsky coefficient [50, 
l5l| . (• • • ) a Clebsch Gordon coefficient and $„im(r) the 
single particle harmonic oscillator wave function. 



$nim(r) ^ R„i{r)Yim{r) 



(Bll) 



with 



Rni{r)=K^A -!-tL('+V2)('ziV-K^)' (B12) 



flho 



*ho 



and 



/ 2 n! aj7^ 



T{n + l + 3/2)' 



(B13) 



with 



^" 



*! «ho 



r(i + 3/2)V2' 



(B5) 



In Eq. (IB4p . the L^- denote the associated Laguerre 
polynomials. Using the generating function of the con- 
fluent hypergeometric function [29t] . 

oo ^(1/2), N 

r(-g)t/H,|,x)=^i^^^^, (B6) 

the interacting wave function V'Joo('"i2) can be ex- 
panded in terms of the non-interacting wave functions 

/ rel.ni / \ 



V'roo(ri2) = EC.V'-oo"(ri2), 



i=0 



where 



a = 



1 



Nf^^^i-q)i^~q)■ 



(B7) 



(B8) 



Inserting the right hand side of Eq. (jB7p into Eq. (|B1 
the one-body density matrix reads 



oo oo 



,(r'i,ri) = 25]5]C,^a 



(B9) 



/ 



i=0 i'=0 

[^SS.(R'i2)Cm>'i2)]*V'ooo(Ri2)V'Ioo'">i2)d'r2. 



Step (ii): To facilitate the integration over r2 in 
Eq. (jB9[) . we expand the product of the non- interacting 



In Eq. (jBlOp . A denotes the total angular momentum 
quantum number to which the two-particle state on the 
left hand side is coupled and II the corresponding pro- 
jection quantum number. For the state of interest, we 
have A = since L = I = 0. Correspondingly, we have 
n = 0. This implies that the sums on the left hand 
side of Eq. (jBlOp reduce to a single term with Clebsch- 
Gordon coefficient (0,0, 0,0|0,0) = 1. For A = B = 0, 
the Clebsch-Gordon coefficient on the right hand side of 
Eq. (jBlOp is only non-zero if I2 ~ h and 1712 ~ —mi, 
which eliminates the sums over I2 and TO2 and yields 
(Zi,mi,Zi,-mi|G,0) = (-l)'i-™i (2^1 + l)-i/2. Using 
these constraints for the quantum numbers, the Talmi- 
Moshinsky bracket on the right hand side of Eq. (IBIOP 
reduces to 1521 



((ni,Zi,n2,;i;0|0,0,z,0;0)) = 



2^ ^-^^ ■ 'ni\n2\ KlNl^ 



4)ii 



-^v/217 + T- 



(B14) 



Energy conservation implies that i is constrained to take 
the values i — ni + n2 + h in Eq. (JB14I) . Applying 
Eq. (jBlOp twice to the integrand of Eq. ((B9l) . with the 
associated restrictions on the quantum numbers, we find 



P(r'i,ri) = 2 



E E 

n'^l'^'m'^n2 nilimin2 



[(2;i + l)(2Z'i + l)]-i/2x 



iCn{+n'2+l[)*Cni+n2+li {-^) ^ '"' ' '"'X 

{{ni,li,n2,li;0\0,0,ni + ^2 +/i,0;0)) x 
{{n[J[,n'^,l[;0\0,0,n[+n'^+l[,0;0)) x 

[*n;iimi(l-'i)]*$„iiimi(l"l)x 

[^n',l[~m[{r2)]*<^n2h-rm{r2)d^r2, (B15) 
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where the sums over i and i' have been ehminated due 
to the energy conservation constraint. 

Step (iii): The integration over r2 only gives non- 
vanishing contributions if n'2 = n2, I'l ~ h and m[ = rrii. 
We thus obtain 



Kr'i,ri)=2 



where 



J2 IIc^\„j*"i'i'"iW)]**"i'i™i('"i)' 

(BI6) 



c\ = y 

n2— 



2Zi + l 



((ni,?i,n2,/i;0|0,0,ni +"2 + ^i,0;0>>x 

((n'l, ?i, 712, /i; 0|0, 0, n'l + nj + ^i, 0; 0)>. (B17) 

The projected one-body density matrix p\^J^(r\^r\)^ 
Eq. (IT^ . can now be calculated readily. In the following 
we consider the case where A = and drop the super- 
script of c ^ for notational convenience. We find 



and 



C, 



(4) 



:(l2/io4 -I- 48^3/10,1 + 17^0,; 



128772 



and, for i > 0, 



Cf =0, 



and 



sp 



^(1) ^ ^00 



1 



c: 



^30 \jV2TrJ 
1 



(2) 



jsp 



N'o \jV2^ 



(1- j^oa)' 



C 



(3) 



-'^00 



^^■0 VjV27r 



3r 



1 - 2j7lo,l 



J (3/^0,2 + 7/i^,i) 



(B24) 

(B25) 
(B26) 

(B27) 



(B28) 



Poo 



{r[,ri) = 2 ^ Cn'^niRn{Q{ri)Rmo{ri), (B18) 



n^ni 



The hn^p are defined in Eq. (IT^ . Using the notation 
introduced in Eq. (jA14p . we find 



where the c„' „j can be interpreted as elements of a 
symmetric coefficient matrix whose eigen values are the 
scaled occupation numbers n^oo/2. The n,yoo/2 are 
shown in Fig. [S] 

In the weakly-interacting regime, we obtain analytic 
expressions for the occupation numbers of the ground 
state by expanding around q — 0. Using Eq. (117^ with 
Te = 0, we rewrite the c„'„j in terms of a; = as(0)/aho 
as opposed to q. Our goal is to obtain the condensate 
fraction of the weakly-interacting two-body ground state 
up to fourth order in x. Extending the analytical pro- 
cedure discussed in Appendix [XJ this requires that we 
calculate cqo up to fourth order in x, cjq up to third or- 
der, and Cjj up to second order. Inspection of Eq. (IB17|) 
shows that the as(0)-dependence of c„' m comes from the 
Cni+n2 and (C„'+„2)* coefficients. We write 

C, « Cf + C^'^x + Cf^x^ + Cf\^ + Cfx^ + O(x^). 

(B19) 

(k) 

The C s needed to evaluate the condensate fraction up 
to order a;"* are 



(2) 
72 



(3) 
72,1 



-2C, 



0'^ - ^ 4^3(1, 1, 1, 5/2, 2, 2, 2, 1/4), (B29) 



5^4(1,1,1,1,5/2,2,2,2,2,1/4) 



4(27r)3/2 



^01 hoih 



02 



«03 



6(27r)3/2 2(27r)3/2 3(27r)3/; 



■j(B30) 



00 00 



7^'2^ = EE 
]=i k=i 



(3) 
72,3 



2l/2-2j--2fcp(^- ^ fc ^ 3^2) 

4jfc(j + fc)7r2 j! fc! 



f. (2) 
Vl72 . 



/27r 



and 



72 



(4) 



0.406786416075. 



(B31) 



(B32) 



(B33) 






(1) 









(3) 



3(27r)3/2 



(B20) 
(B21) 

(B22) 



(/io,3 + 3/10,2/10,1 +2/ii^,i), (B23) 



In Eqs. (|B29|) and (|B30[) . qFp denotes the generalized 

(2) 
hypergeometric function. The numerical values of 72 , 

72 1 ' 72 2 and 72 3 are listed in Table IIIII 

The approach discussed above can be extended to ac- 
count for the effective range dependence of the conden- 
sate fraction, yielding the result discussed in the last 
paragraph of Sec. IIIII 
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